#!/bin/python
#-*- coding: utf-8 -*-
import struct
import numpy as np
import os
import gdal
import osr


def parse(rfe_path):

	if rfe_path.endswith(".gz"):
		os.system("gunzip -k %s"%rfe_path)

	rfe_file = open(rfe_path.replace(".gz",""),'rb')

	data = rfe_file.read()

	# in south asia, the wdith and height of image are 401 and 301
	raster_array = np.ones((301,401)).astype(np.float)

	for i in range(0,301):
		for j in range(0,401):
			raster_array[300-i][j]=struct.unpack('!f',data[(i*401+j)*4:(i*401+j+1)*4])[0]

	rfe_file.close()
	os.remove(rfe_path.replace(".gz",""))

	return __add_geo_info__(raster_array)

def __add_geo_info__(raster_array):
	out_raster = gdal.GetDriverByName("MEM").Create("",raster_array.shape[1],raster_array.shape[0],1,gdal.GDT_Float32)

	projection = osr.SpatialReference()
	projection.SetWellKnownGeogCS('WGS84')
	geotrans = (70,0.1,0.0,35.00,0.0,-0.1)

	out_raster.SetGeoTransform(geotrans)
	out_raster.SetProjection(projection.ExportToWkt())
	out_raster.GetRasterBand(1).SetNoDataValue(-999.0)

	out_raster.GetRasterBand(1).WriteArray(raster_array)
	return out_raster

if __name__ == '__main__':
	parse("/mnt/fire/Metero/rfe/cpc_rfe_v2.0_sa_dly.bin.20130102.gz")
